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Abstract 

Background: The complement protein C5a acts by primarily binding and activating the G-protein coupled C5a 
receptor C5aR (CD88), and is implicated in many inflammatory diseases. The cyclic hexapeptide PMX53 (sequence 
Ace-Phe-[Orn-Pro-dCha-Trp-Arg]) is a full C5aR antagonist of nanomolar potency, and is widely used to study C5aR 
function in disease. 

Results: We construct for the first time molecular models for the C5aR:PMX53 complex without the a priori use 
of experimental constraints, via a computational framework of molecular dynamics (MD) simulations, docking, 
conformational clustering and free energy filtering. The models agree with experimental data, and are used to 
propose important intermolecular interactions contributing to binding, and to develop a hypothesis for the 
mechanism of PMX53 antagonism. 

Conclusion: This work forms the basis for the design of improved C5aR antagonists, as well as for atomic^detail 
mechanistic studies of complement activation and function. Our computational framework can be widely used 
to develop GPCR-ligand structural models in membrane environments, peptidomimetics and other chemical 
compounds with potential clinical use. 

Keywords: Class A GPCR, C5aR, C5a, Complement system, Molecular dynamics. Docking, Implicit solvent. 
Membrane protein 



Background 

C5aR is the membrane-bound receptor for the comple- 
ment system protein C5a and the target of inhibition 
against inflammatory diseases. Here, we present the 
development of a structural model for the complex of 
membrane-embedded C5aR and its antagonist peptide 
PMX53. We accomplish this by an innovative and com- 
prehensive computational framework that combines con- 
formational sampling for both receptor and ligand with 
docking, and evaluates a large number (-300,000) of 
docked conformations by structural and free energy cri- 
teria. Even though we do not impose any experimental 
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restraints, the resulting structural models are consistent 
with available experimental data. 

The complement system is a major and essential com- 
ponent of the innate immune response. It can be acti- 
vated following infection or injury through four distinct 
pathways, which lead to opsonisation of pathogens, cell 
lysis, and the production of potent pro-inflammatory 
peptides. The complement protein C5a is generated fol- 
lowing cleavage of the 5* component of complement 
(C5), and is one of the most potent inflammatory media- 
tors in humans [1]. Given its potent inflammatory activ- 
ity, prolonged or inappropriate activation of complement 
can generate unwanted C5a, which is implicated in many 
inflammatory diseases [2,3]. 

C5a induces the majority of its known effects primarily 
through a G-protein coupled receptor, termed C5aR. It 
can also bind to a second receptor called C5L2, however 
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this receptor does not couple G-proteins, and thus has 
unclear and controversial functions [4]. The C-terminal 
ten residues of C5a have been shown to be critical for 
C5aR activation, and can induce full activation (efficacy) 
of C5aR, albeit at lower potencies than intact C5a. This 
C-terminus is proposed to act at the second extracellular 
loop of C5aR [5]. 

Given the proposed involvement of C5a in classical in- 
flammatory diseases such as rheumatoid arthritis, [6] and 
a widening of the pathogenic roles of C5a to traditionally, 
non-inflammatory diseases such as cancer [7], and brain 
diseases [8], there has been a burgeoning interest in devel- 
oping inhibitors of C5a for clinical use [9]. One such 
method has been to block C5a-C5aR interaction through 
the development of selective C5aR antagonists. In 1999 
and 2000, a series of cyclic peptide full antagonists were 
described which were designed around the C-terminal res- 
idues of C5a [10,11]. One of these compounds, termed 
3D53 (later named PMX53), has been the most exten- 
sively studied, and is widely used to study C5aR function 
in disease [3]. It is a cyclic hexapeptide (structure: Ace- 
Phe-[Orn-Pro-dCha-Trp-Arg]), which has low nanomolar 
potency, and efficacy in a wide range of species [12], mak- 
ing it an ideal research tool [3]. One potential drawback 
for the clinical use of this compound, however, is its pep- 
tidic nature and low oral bioavailability [13], which has 
limited its commercial development [9]. 

The recent determination of several GPCR crystal struc- 
tures has contributed significantly to our understanding of 
GPCR structural organization, ligand specificity, and acti- 
vation [14-19]. MD simulations have provided consider- 
able insight on the conformations of inactive states and 
the mechanism of receptor activation (Ref [20] and refer- 
ences therein). Nevertheless, the generation of reliable 
structural models for GPCR proteins of unknown struc- 
ture and for their ligand complexes remains a challenging 
task [21]. 

In the present study we develop molecular models of 
PMX53 bound to C5aR in a model membrane. Our goal 
is to use the resulting models to develop peptidomimetics 
[22] and chemical compounds, which may overcome 
some of the clinical hurdles that were encountered with 
PMX53 [3]. We start from a well-tested earlier structural 
model for free C5aR [23-27] and an NMR structure for 
PMX53 [28], and generate a large number of representa- 
tive conformations for free C5aR and the PMX53 ligand 
by multi-ns MD simulations. Docking of the obtained 
conformations yields a large number of structural models 
for the complex (-300,000 structures); subsequent filter- 
ing by structural and energetic criteria identifies a man- 
ageable number of candidate complexes, which are 
further tested by extensive MD simulations and free 
energy calculations. In contrast to previous computational 
studies of the C5aR complexes with C5 and PMX53, we 



represent membrane effects on the protein and ligand in- 
teractions at the MD simulation and assessment stage via 
high-quality implicit-membrane models [21,29-31]. Such 
models offer a promising alternative to explicit-membrane 
treatments, which are more accurate but are impractical 
for structure-prediction calculations [32,33]. Furthermore, 
the implicit-membrane representation enables the rapid 
determination of binding free energies for several hundred 
thousand candidate structures. Even though no experi- 
mental information has been incorporated a priori in the 
docking, the most promising complexes are consistent 
with available experimental data, reflecting the accuracy 
and potential of the employed methodology. 

The obtained models for the C5aR:PMX53 complex 
can serve as the basis for knowledge-based discovery of 
C5aR antagonists with improved properties compared to 
PMX53, as well as for basic mechanistic studies of com- 
plement activation and function at molecular detail and 
atomic resolution. Furthermore, the described combin- 
ation of implicit-membrane MD simulations, docking and 
free energy calculations is a promising framework for 
the generation and assessment of structural models for 
GPCR-ligand complexes. 

Methods 

Description of simulation systems 
C5aR receptor 

The human GPCR receptor C5aR consists of 350 amino 
acids, and has the typical GPCR topology, with an extra- 
cellular N-terminal fragment, seven trans-membrane 
(TM) helices interconnected by extracellular (EC) and 
intracellular (IC) loops, and an intracellular C-terminal 
fragment [26]. Nikiforovich et al. has constructed struc- 
tural models for free C5aR [24] as well as its complex 
with C5a [25,26]. In the MD simulations we use as a 
starting point for C5aR the structural model of Nikiforo- 
vich [24-26]. The seven transmembrane helices in the 
Nikiforovich model are defined as: 38-63 (HI), 71-98 
(H2), 107-138 (H3), 150-172 (H4), 199-224 (H5), 236- 
267 (H6), and 281-300 (H7); similarly, the three extra- 
cellular loops are defined as 99-106 (ECl), 173-198 
(EC2), and 268-280 (EC3) [26]. The simulation system 
omits the first seven amino acids, which do not contrib- 
ute to C5a binding and are not expected to affect bind- 
ing of PMX53. It also omits the intracellular C-terminal 
region 307-350, which is very remote from the insertion 
point of the C5a C-terminal end, and the putative ligand 
binding site. 

PMX53 ligand 

The hexapeptide PMX53 (Figure 1) has the sequence Ace- 
Phe-[Orn-Pro-dCha-Trp-Arg]; Ace denotes the blocking 
group CH3-CO at the N-terminal end, Orn ornithine, dCha 
d-cyclohexyl-alanine, and the brackets denote cyclization of 
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Figure 1 Structure of the cyclic hexapeptide ligand PiVlXSB, with amino acid sequence Ace-Phe-[Orn-Pro-dCha-Trp-Arg]. Panel A: 

two-dimensional Clnemical structure PMX53. The N-terminal end is blocked by the group CH3-C0 (Ace); Orn denotes ornithine and dCha 
d-cyclohexyl-alanine. The brackets denote cyclization via a covalent bond between the Orn side-chain and the Arg6 carbonyl group. Panel B: 
three-dimensional representation of the first conformer of the NMR ensemble of structures of PMX53. Atoms are colored according to element 
type. Panel C: three-dimensional representation of the ensemble of the NMR structures of PMX53. The color of residues transitions from blue to 
red between the N- and C- termini. Hydrogens are omitted for clarity. 



the mainchain via a covalent bond between the Orn side- 
chain and the Arg6 carbonyl group. Figure lA shows 
the chemical structure of PMX53, and Figure 1B,C shows 
three-dimensional representations of the NMR structure 
of PMX53 [28]. 

Construction of structural models for the C5aR:PMX53 
complex 

In order to achieve a systematic construction and evalu- 
ation of plausible structural models for the complex, we 
employed a range of methods, including MD simulations, 
docking, energy minimizations, and binding-affinity calcu- 
lations. Our computational framework consisted of the 
following steps: (A) generation of a large collection of rep- 
resentative PMX53 and C5aR structures via long MD sim- 
ulations of the isolated ligand and receptor; (B) clustering 
of the simulation trajectories and determination of high- 
probability conformations; (C) generation of a large num- 
ber of potential structural models for the complex, 
via docking of conformations from the most populated 
PMX53 and C5aR clusters; (D) filtering of the models with 
structural and energetic criteria; (E) assessment of the most 
promising models by MD simulations and binding free en- 
ergy calculations. Each step is detailed below. 



the embedding membrane. Hence, it is important to gen- 
erate alternative structural models of PMX53, which are 
likely to be stabilized by the less polar environment of 
the binding site. As an approximation, we generated such 
conformations by MD simulations of free PMX53 in a 
membrane environment, represented by the implicit- 
membrane molecular-volume Generalized Born (GBMV) 
model [31]. Standard amino acids were described by 
the CHARMM all-atom topology and energy function 
[34], using the CMAP correction for all L-amino 
acids [35]. Topology and parameters for non-standard 
groups (ornixthine, cyclohexylalanine and the Orn2-Arg6 
cyclization segment of PMX53) were derived, respect- 
ively, from CHARMM22 definitions for lysine, alanine/ 
cyclohexane and the peptide group. 

To represent conditions of variable polarity, we con- 
ducted four simulations, in which we restrained the 
PMX53 center of mass at 0 A, 10 A, 15 A and 20 A, re- 
spectively, from the membrane center. The total simula- 
tion length was 180 ns (45 ns per position). We also 
conducted an additional 15-ns run in aqueous solution, 
represented by the corresponding aqueous GBMV model 
[30]. All simulations were conducted with the CHARMM 
program, version c35b3 [36] . 



Generation of PMX53 conformations 

Competition binding studies with linear and cyclic pep- 
tide antagonists suggest that the binding site of PMX53 
is in the transmembrane region of C5aR, near or at the 
location of the binding site of the C5a C-terminal end 
[5]. NMR experiments [28] have shown that the domin- 
ant conformation of free PMX53 in deuterated DMSO 
(DMSO-d6) has residues 1-2 in a random-coil state, and 
segment 3-6 in a type-II |3-turn. PMX53 may assume a 
different conformation in the complex with C5aR, due 
to the influence of the surrounding C5aR receptor and 



Classification of PMX53 conformations 

The dominant simulation conformation of the PMX53 
mainchain was in agreement with the NMR structure 
[28], with residues 1-2 in random coil and segment 3-6 
in a type-II P-turn. We partitioned the conformations 
into groups, based on the following nine side-chain dihe- 
dral angles: (i) IN-lCa-lCp-lCy, (ii) ICa-lCp-lCy-lCsi, 
(iii) 4N-4C„-4Cp-4C^ (iv) 4C„-4Cp-4CY-4C5i, (v) 5N- 
5Ca-5Cp-5Cy, (vi) 5Cc, -5C|3-5Cy-5Csi, (vii) 6N-6Ca-6Cp- 
6Cy, (viii) 6Ca-6C|3-6Cy-6C5, (ix) 6C|3-6Cy-6C6-6Ni;. We 
employed the clustering algorithm of the CHARMM 
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program [37,38], a clustering radius of 45° and a "max- 
imum error" radius (reflecting the maximum difference 
among the centers of distinct clusters) of 10°. We ob- 
tained a total of 775 clusters, with 51 clusters containing 
ca. 30% of the conformations. We employed represen- 
tative structures from all 51 clusters in the docking 
calculations. The clustering focused on side-chain di- 
hedral angles as the side-chains are significantly more 
flexible than the backbone. For the 51 clusters, the average 
root-mean-square difference (RMSD) of the backbone/ 
side-chain heavy atoms with regard to the average 
structure of the 51 clusters is equal to 0.54 ± 0.15 A/ 
2.46 ± 0.52 A. A significantly higher mobility of the 
side-chains compared to the backbone is also present 
in the 10 NMR-derived structures (see Figure IC), ac- 
cording to which the average RMSD of the backbone/ 
side-chain heavy atoms with regard to the average NMR 
structure of the 10 conformations is equal to 0.23 ± 
0.07 A/1.85 ±0.22 A [28]. 

Generation of C5aR conformations 

The initial C5aR conformation was taken from the Niki- 
forovich model for the C5aR:C5a complex [25]. Details 
on the construction of the structural models for C5aR 
and the C5aR:C5a complex can be found in previous 
work [25,39]. 

Starting from this structure, we generated conforma- 
tions of C5aR via MD simulations with the implicit- 
membrane switching-function generalized Born (GBSW) 
module [29]. Two parameters of this model are the total 
thickness of the low-dielectric membrane slab {Tynemb) 
and the half-length {Imsw) of the membrane/water inter- 
face buffer region. Based on UNIPROT [40] definitions, 
we estimated the total thickness of the intra-membrane 
C5aR helical region at ~ 35.5 A. At the same time, the 
OPM database [41] yields a hydrophobic thickness of 
32.2 ± 1.2 A for rhodopsin (PDB code IGZM), the pro- 
tein on which the C5aR model is based. We combined 
the OPM and UNIPROT predictions, by setting T^^mZ, 
to 36 A and l^sw to 2.5 A. With these values, the com- 
bination Tynemb " 2* l^sw IS 31 A, 1.6. correspouds to the 
lower bound of the OPM estimate. Test simulations 
indicated a stable behavior (RMSD of ~3.0 A from 
the initial conformation), whereas somewhat different 
Tmemb values resulted in larger RMSD from the initial 
conformation (~3.5 A). A similar membrane thickness 
(35 A) and half-length (2.5 A) were used in recent 
structural modeling of class A GPCRs, which employed 
the same implicit-membrane model [21]. The surface 
tension coefficient (y) was set to 0.04 kcal/(mol A^); 
other parameters were set to default values of the GBSW 
model [21]. 

To enlarge the ensemble of generated C5aR conforma- 
tions, we conducted simulations with four protocols: In 



protocol (i) we removed the C5a ligand and simulated 
the C5aR protein at an elevated temperature (400 K), 
with harmonic restraints on the heavy backbone atoms 
of the entire protein, and the side-chains unrestrained. 
In this way, the C5aR binding pocket retained a simi- 
lar volume as in the C5aR:C5a complex (Nikiforovich 
model) [25], while the simulation eliminated any bias in 
the initial side-chain conformations, which might arise 
due to interactions with the C5a ligand. In protocols ii-iii 
we simulated the complex between C5aR and the C5a 
fragment 60-74, which binds at the intra-membrane 
C5aR region. The C5a fragment was either restrained 
near its initial shape via the bestfit option of CHARMM 
(protocol ii), or was left unrestrained (protocol iii). Fi- 
nally, in protocol iv we simulated the entire C5aR:C5a 
complex without restraints. The simulation temperature 
of protocols ii-iv was 300 K. Prior to each production 
run, we subjected each system to energy-minimization, 
heating, and 1.6-ns equilibration. The length of the pro- 
duction run was 11 ns for protocol iii, and 5.5 ns for all 
other protocols. 

Standard amino acids were described by the CHARMM 
all-atom topology and energy function [34] including a 
CMAP correction of the backbone torsional angle energet- 
ics [35]. A 16-A cutoff distance was used for non-bonded 
interactions. The lengths of covalent bonds containing 
hydrogen atoms were constrained by the SHAKE algo- 
rithm [42], and the equations of motion were solved with 
an integration time step of 2.0 fs. The system was in con- 
tact with a Langevin heat bath at 300 K; a friction coeffi- 
cient of 5 ps" was used for heavy atoms. 

Classification of CSaR conformations 

We characterized the C5aR conformations by the shape 
of their intra-membrane binding pocket, since structural 
differences in other protein regions should not be as 
relevant to PMX53 binding. We described the shape of 
the binding pocket by a novel methodology, which filled 
the binding cavity (after removal of C5a) with particles 
having the approximate diameter of a water molecule 
and taken from a water box; these particles created a 
grid of points, with an inter-point spacing equal to the 
water diameter. The underlying premise was that by en- 
suring that the particles were always in the same pos- 
ition relative to C5aR, it would be possible to identify 
the cavity regions that were changing structurally, by 
simply keeping track of which particles fit in the cavity. 
To create an ensemble of cavity shapes, we first super- 
imposed each MD snapshot onto the initial structure 
(the Nikiforovich model), based on the Ca-atoms of the 
C5aR trans -membrane regions; we then overlayed a 
50 A X 50 A X 60 A explicit water box on each snapshot, 
and deleted all water molecules whose oxygen atom was 
within 2.4 A of any C5aR heavy atom. Even though C5a 
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was not taken into account when computing the cavity 
shapes, water molecules within 4 A of the coordinates of 
the 15 C-terminal residues (60-74) of C5a were used to 
define the top of the cavity; water molecules further than 
4 A from these residues were removed. 

We created lists identifying the remaining molecules, 
and used them to generate binary water fingerprints 
with length equal to the number of unique water mole- 
cules found in the cavities of the MD snapshots. A value 
of 1 (0) was associated with present (absent) water mole- 
cules. We quantified the similarity between two water 
fingerprints A and B via distance matrices based on the 
Jaccard binary distance measure [43]. 

r _ '-'10 ~'~ \ 

Jab - r 4-r 4-r ^ ' 

In Eq. (1), Cj2 is the number of common water mole- 
cules in both fingerprints (cavities), and Cjo. Coi are, re- 
spectively, the numbers of water molecules only in 
fingerprint A or B. The Jaccard distance varies from 0 
(fingerprints with the same water molecules) to 1 (no 
water molecules in common). We then used hierarchical 
clustering [44] to classify the MD snapshots into fam- 
ilies, based on water fingerprint similarity, and visualized 
the clustering via dendrogram trees. The dendrogram 
trees were cut at an arbitrary binary similarity value of 
0.3, to keep the number of selected structures manage- 
able. With this value, two structures not in the same 
cluster differed by more than 50 water molecules. From 
each cluster, we selected as most representative the cav- 
ity structure with the largest volume (largest number of 
present water molecules). Ultimately, we chose 150 - 
250 structures from each protocoU (a total of 785 struc- 
tures) for the docking studies. These calculations were 
performed with scripts written in the R statistical lan- 
guage [45] using the Bio3D [46] library. The above 
methodology was not employed to the high-temperature 
simulation (protocol i), which was restrained to maintain 
the volume of the binding cavity; from that trajectory, 
we simply extracted structures at intervals of 40 ps. 

Docking 

Using representative C5aR and PMX53 snapshots from 
the structural classification described above, we generated 
possible structures of the C5aR:PMX53 complex with the 
program DOCK6 [47]. For each C5aR snapshot, we cre- 
ated the corresponding molecular surface via DMS [48], 
and represented potential binding sites by spheres via the 
SPHGEN utility. SPHGEN identifies spheres tangent to 
the molecular surface, and performs clustering to elimin- 
ate redundant spheres, which represent surface cavities 
and serve as potential locations of atoms during ligand 
orientation. For each C5aR snapshot combination, we only 



retained spheres within 7 A of the coordinates for the 15 
C-terminal residues of C5a. For protocols ii-iv, the coordi- 
nates for C5a were taken from the specific snapshot; for 
the elevated temperature run, the C5a coordinates were 
taken from the Nikiforovich model. A receptor-ligand 
clash was defined as >50% atom-atom overlap of two 
atomic radii; we excluded docking poses with more than 
five such clashes. All orientations passing this filter were 
saved as input for the CHARMM-based scoring proced- 
ure, since the docking program was not used to rank the 
C5aR:PMX53 poses, due to the complexities of GPCR 
docking. 

Altogether, 51 PMX53 conformations, corresponding 
to the most populated 51 clusters from the PMX53 tra- 
jectory analysis, were docked to the representative struc- 
tures from each of the four C5aR MD simulations, 
resulting in 306,497 structural models for the complex. 
A rigid docking procedure was used; nevertheless, the 
use of flexible templates of the ligand and the receptor 
incorporate flexibility in the docking. 

Filtering of clocked structures 

We subjected all 306,497 docking conformations to 
100 energy steepest-descent minimization steps in a 
uniform dielectric medium (e = 4), followed by 50 steepest- 
descent minimization steps in the GBSW implicit- 
membrane environment. At the end of minimization 
we filtered the conformations by a combination of energy- 
based and structure-based criteria. A flowchart is pre- 
sented in Figure 2. 

The energy-based criterion was applied in two steps: 
First, we computed the binding free energies of all 
306,497 conformations (at the end of minimization) in 
the Molecular Mechanics Generalized Born Surface Area 
Approximation (MM-GBSA) [49-51], via the following 
equation 

AG = Gpl-Gp-Gl (2) 

In Eq. (2), Gx is the total free energy of molecule X 
(complex PL, free protein P, or free ligand L). The pro- 
tein and ligand conformations were assumed identical 
in the complex and free molecules. With this assump- 
tion, any bonded-energy contributions to AG cancel in 
Eq. (2). Even though this assumption is not strictly accur- 
ate, it is probably sufficient for the identification and 
rejection of conformations with weak association free 
energies. Protein, ligand and complex entropic contribu- 
tions are ignored; they are expected to approximately 
cancel when comparing relative affinities of different 
binding modes. The solvation components of the com- 
plex and free protein free energies were computed in the 
inhomogeneous membrane/water environment; the solv- 
ation component of the unbound ligand free energy was 
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MD simulations of free C5aR & 
clustering 785 structures 



MD simulations of free PMX53 & 
clustering -^51 structures 



Docking calculations -> 306947 binding 
models 



i 



Energy minimization and evaluation of binding free energies (GBSA 
approximation) ^ 4676 structures with free energies between -67 and -40 kcal/mol. 



i 



Re-evaluation of binding free energies (PBSA approximation) & structural clustering (based on 
ligand main-chain coordinates) -> 35 distinct binding modes with lowest binding free energy 



i 



Assessment of binding modes by MD simulations and PBSA binding free-energy 
calculations: Comparison with experimental data 



Figure 2 Flowchart of the computational framework used to generate and assess structural models of the C5aR:PMX53 complex. 



computed in aqueous solution (modeled by the impli- 
cit GBSW model) [29]. We thus identified 4676 struc- 
tures (resulting from the minimization of the 306,497 
docking conformations), whose binding free energies were 
within 27 kcal/mol of the binding free energy minimum 
(-67 kcal/mol). 

In the second filtering step with the energy-based 
criterion, we calculated the binding free energies of these 
4676 structures in the Molecular Mechanics Poisson- 
Boltzmann Surface Area Approximation (MM-PBSA) 
[49]. In these calculations, we inserted the complex and 
protein in a membrane slab with a thickness of 31 A 
and a dielectric constant of 1, surrounded by water with 
a dielectric of 80. The protein dielectric constant was 
set to 2. The free ligand was placed in pure water. The 
MM-PBSA calculations were performed with the Poisson- 
Boltzmann solver of the CHARMM program (PBEQ 
module). We used 150 grid points in each direction and a 
grid-spacing of 0.5 A. The MM-PBSA calculation yielded 
a more expanded free energy range (between -90 kcal/ 
mol and -20 kcal/mol), and changed the relative affinity 
for some structures. The correlation between MM-GBSA 
and MM-PBSA was analyzed via a linear least squares fit; 
the slope and the standard error, respectively, were equal 
to 1.03 and 0.03. 

In the next step, we used a structure-based criterion 
in conjunction with the MM-PBSA binding free ener- 
gies, to identify structurally distinct, high-affinity binding 
modes. Clustering analysis of the PMX53 mainchain coor- 
dinates (N, CA and C atoms) with a 5-A clustering radius 
identified 35 clusters, each reflecting a distinct ligand- 
binding mode. The analysis was performed with the 
program WORDOM [52]. From each cluster we extracted 



the conformation of lowest MM-PBSA binding free en- 
ergy, and used it as a starting point in MD simulations. 

Evaluation of structural models for the C5aR:PMX53 
complex by MD simulations 

Prior to the MD runs, we relaxed each structure via 
1000 steps of steepest-descent energy minimization. 
Subsequently, we heated the complexes to a temperature 
of 300 K via four 100-ps runs with respective tempera- 
tures of 75 K, 150 K, 225 K and 300 K; during the heat- 
ing stage we restrained all heavy backbone atoms to 
their initial positions by a harmonic force with strength 
of 5 kcal/mol/A^. We then conducted three 400-ps 
equilibration runs at 300 K, in which we gradually re- 
duced the harmonic-restraint force constant to 1 kcal/ 
mol/A for intra-membrane residues, and to zero for 
residues in extra-cellular (EC) and intra-cellular (IC) 
loops. In the subsequent production runs, we removed 
all protein and ligand restraints. Parameters and simula- 
tion details were the same as for the free ligand and 
C5aR simulations, described above. The duration of each 
production run was 7 ns, with only the last 4 ns employed 
in the analysis. 

In order to assess the stabilities of the simulated com- 
plexes, we computed the corresponding MM-PBSA as- 
sociation free energies. For each run, the results were 
averaged over 700 snapshots extracted at 10-ps intervals. 
The best complex had a free energy of -175 kcal/mol. 
Seven complexes within 20 kcal/mol of the lowest value 
(-175 kcal/mol) of less than -155 kcal/mol (20 kcal/mol 
greater than the average of the most promising complex) 
were simulated for an additional 13 ns (a total of 20 ns); 
the last 17 ns were employed in the analysis. 
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Residue pairwise interaction free energies 

We analyzed the intermolecular interaction free energies 
between PMX53 and C5aR residue pairs of the final 
simulation using the following equation: 

i^R j^R' i^R jeR' i^R^R' 

» ^polar A ^non polar 

^^RR' ^^RR' 

(3) 

The first and second group of terms on the right-hand 
side of Eq. (3) describe, respectively, polar and non-polar 
interactions between a C5aR residue R and a PMX53 
residue R'. The polar component contains a Coulombic 
energy term, and a free energy GB contribution, model- 
ing the interaction between residue R and the solvent 
polarization potential induced by R'. Similarly, the non- 
polar component contains a van der Waals interaction 
between R, R' and a surface free energy term, expressing 
cavity contributions and nonpolar interactions with the 
surrounding solvent. 

The non-polar and polar solvation terms were calculated 
using the implicit membrane GBSW. The generalized- 
Born energies and the atomic accessible-surface areas (ASj) 
entering in Eq. (3) depend on the location of R and R' in 
the complex. To compute the polar (GB) interaction free 
energy term for a particular residue pair (R, R') in Eq. (3), 
we set to zero the charges of atoms not in residues R and 
R'. The surface term contains the difference in solvent ac- 
cessible surface areas of residues R and R' in the complex 
and unbound states; the interaction free energy non-polar 
term represents the creation of a cavity in its surrounding 
medium (membrane/water) to accommodate biomolecules 
and switching-on dispersion interactions between biomole- 
cules and the surrounding medium while all atomic 
charges are set to zero [51]. A similar methodology for the 
analysis of interacting residues has been used for the eluci- 
dation of species-specificity of complement protein C3 
[53] and the design of transgenic (mouse/human) C3 with 
putative affinity for compstatin [54] . 

Results and Discussion 

Earlier models of free C5aR, and the complex of C5aR 
with its natural activator, peptide C5a 

The C5aR protein has the typical GPCR topology, with an 
extracellular N-terminal fragment, seven trans-membrane 
(TM) helices interconnected by extracellular (EC) and 
intracellular (IC) loops, and an intracellular C-terminal 
fragment [55]. Nikiforovich et al. has constructed struc- 
tural models for free C5aR [24], as well as its complex 
with C5a [25,26]. The structure of the C5aR TM re- 
gion in these models was based on the corresponding 
TM region of dark-adapted rhodopsin, and the loops 



were constructed with a de novo structure prediction 
method [23]. The model has been used to interpret ex- 
perimental results on C5a binding to C5aR [25], and to 
model the conformational changes occurring during re- 
ceptor activation in the TM region [56] and EC loops [27] 
of C5aR. The loop prediction method was also applied 
with success to construct structural models of extracellu- 
lar loops for other GPCR receptors [57]. 

Numerous systematic experimental studies have led to 
a two-site model of C5aR activation (see [2,3] and refer- 
ences therein). The primary affinity site involves contacts 
between acidic residues in the N-terminal end of C5aR 
and basic residues in the core of C5a; interactions in this 
site contribute to binding strength, but not to C5aR acti- 
vation. The second site is formed by residues in the TM 
domain and the EC loops and interacts with the C5a ter- 
minal fragment 60-74. This site contributes to the acti- 
vation of C5aR by C5a and other peptide agonists. 

Experiments with agonist and antagonist peptide mimics 
and receptors mutated at possible interaction sites have 
provided information on the ligand binding site of C5aR. 
Competition-binding experiments suggested that PMX53 
binds at or near the same TM location as the C-terminal 
moiety of C5a [5]. The substitution Ilell6Ala correlated 
with the emergence of weak agonist activity in a PMX53 
derivative with a bulkier side-chain (benzothiazolalanine) 
at position 5 [5], and converted a linear peptide mimic 
from antagonist to agonist [58]. These results suggested 
that residue 5 of PMX53 (and related peptides) binds in 
the vicinity of Ilell6. The C5aR substitution Asp282Ala 
(in the ECS loop) caused a 10-fold reduction in PMX53 
affinity, suggesting the formation of an electrostatic inter- 
action between Arg6 of PMX53 and Asp282. The substi- 
tution Argl75Asp converted PMX53 to a weak agonist, 
indicating that Argl75 might play a role in the discrimin- 
ation between agonists and antagonists; additional experi- 
ments suggested a possible interaction between Argl75 
and the C-terminal carboxylate of hexapeptide ligands [5]. 
Mutations Arg206Ala and Glul99Lys also affected the 
receptor activation by hexapeptide ligands; modeling in 
[5] argued that a possible interaction Arg206-Glul99 may 
stabilize the position of helix H5, and disruption of this 
interaction may contribute to activation. 

Taking into account these experimental results, 
Higginbottom et al. constructed two models for the 
PMX53:C5aR complex [5]. The C5aR conformation was 
based on homology with dark-adapted rhodopsin. Loop 
EC2 was kept in a "closed" position (near the membrane), 
via a Cysl09-Cysl88 disulfide bond. Various docking con- 
formations of PMX53 inside the C5aR binding pocket 
were tried, under the condition that they satisfied two 
ligand-receptor interactions: (i) The center of the Trp5 in- 
dole ring was restrained within a distance of 5 A from the 
Ca atom of Ilell6, and (ii) the side-chain of Arg6 formed 
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a salt-bridge with the side-chain of Asp282. In the result- 
ing model, the PMX53 Phel side-chain was positioned be- 
tween helices 6 and 7, in the vicinity of residues Ile263 
(H6) and Leu207, Phe211 (H7); the nonpolar dCha4 side- 
chain was directed toward helices 3 and 4, near Pro 11 3, 
Serll4, Argl97 and the Glul99-Arg206 salt bridge; Trp5 
was placed near the hydrophobic residues Pro 11 3, Ilell6 
and Val286; finally, Arg6 was proximal to Aspl91, Tyr258 
and Gln259, in addition to its salt-bridge with Asp282. 
The final model was presented in Figure five(B) of [5]. 
Higginbottom et al. [5] argued that the activation of C5aR 
may be prompted by the insertion of the Arg6 side-chain 
at a buried position in the ligand binding pocket, past 
Ilell6. On the other hand, the Ar6-Asp282 interaction in 
the modeled C5aR:PMX53 complex prevents the deep in- 
sertion of the Arg6 guanidinium group into the pocket, 
enabling the antagonist activity of PMX53. A second 
model resulted from docking attempts of PMX53 onto 
C5aR with the same two constraints (Trp5-Ilell6, Arg6- 
Asp282), but with the disulfide bridge Cysl09-Cysl88 
absent and the C5aR EC2 loop in an open conformation. 
In the resulting model, (Figure five(C) of Higginbottom 
et al.) [5], the ligand was positioned with the Phel side- 
chain at a similar position (between helices H6 and H7), 
but with the Arg6 and dCha4 interactions swapped. We 
note that the covalent disulfide bridge Cysl09-Cysl88 
should also be present in the open EC2 conformation. 
As we show below, PMX53 is inserted in the second rec- 
ognition region, near the modeled location of the C5a 
C-terminal segment 69-74. 

Ranking of the most promising C5aR:PIVIX53 complexes 

A detailed description of our modeling procedure is pre- 
sented in the Methods, and a summarizing flowchart is 
shown in Figure 2. We first generated a large number of 
representative conformations for the ligand and the recep- 
tor by MD simulations, and grouped them into conform- 
ational clusters. We then performed docking calculations 
of protein and ligand conformations, without any geo- 
metrical restraints. Using a filtering calculation, consisting 
of conformational clustering, energy minimization and 
binding free energy evaluations, we identified 35 promis- 
ing binding modes covering the entire TM cavity of C5aR. 
The MM-PBSA binding-free energies of these modes 
ranged between -90 and -65 kcal/mol (Additional file 1: 
Figure SI). 

We assessed further the stabilities of all 35 structural 
modes by implicit-membrane MD simulations of the 
C5aR:PMX53 complex. In all cases, the initial conforma- 
tions were well preserved; the average RMSD of all TM 
backbone protein residues was 3.3 ± 0.2 A relative to the 
structure of Nikiforovich et al. [25] MM-PBSA analysis 
of the simulation trajectories (up to 7 ns) showed that 
seven structural modes (22, 6, 3, 30, 1, 4, and 10) had 



significantly lower binding free energies. The simulations 
of these models were extended to 20 ns; this length was 
sufficient for convergence of the MM-PBSA association 
free energies (not shown). 

The free energy values (averaged over the last 17 ns 
for the 20-ns trajectories, or last 4 ns for the 7-ns trajec- 
tories) are plotted in Additional file 1: Figure SI as red 
diamonds. The averaged-over-simulation MM-PBSA affin- 
ities are lower than the initial values of the same models. 
This is mainly due to the vdW component, reflecting the 
improvement of protein-ligand contacts with respect to 
the initial docked structures. 

Structural mode 22 has the lowest association free en- 
ergy (-175 kcal/mol). Interestingly, this mode agrees with 
experimental data and with features of an earlier mode 
proposed by Higginbottom et al. [5] The next best modes 
(6 and 3) have association free energies of -166 kcal/ 
mol, ~10 kcal/mol higher relative to mode 22. Mode 
30 has a comparable association free energy with 6 and 3 
(-165 kcal/mol), but involves an entirely different ligand 
orientation. All four modes are described below. 

Structural analysis of mode 22 and comparison with 
experimental data 

The conformation of mode 22 is displayed in Figure 3. 
PMX53 is surrounded by helix 3, helix 4, the extracellular 
domain EC2, helix 5, helix 6, the extracellular domain 
EC3, and helix 7. The N-terminal end of the ligand points 
toward the membrane interior, with the center-of-mass 
(COM) of the Ac-Phel mainchain at a distance of ~ 5 A 
from the membrane center. The C-terminal end of the lig- 
and is directed toward the membrane-water interface; the 
COM of the Arg6 mainchain lies at the interface, at a dis- 
tance of ca. 14 A from the membrane center (Figure 3A). 

A close-up view of the complex conformation in the 
vicinity of PMX53 is shown in Figure 3 (panels B and C). 
Interaction free energy components of all protein-ligand 
residue pairs are displayed in Figure 4, and contact maps 
are shown in Figure 5. Furthermore, a list of PMX53- 
C5aR residue pairs with total intermolecular residue-pair 
energies of at least 1.8 kcal/mol, as evaluated by Eq. (3), is 
included in Table 1. 

Ac-Phel is inserted into a binding pocket formed by 
helices H3, H4, H5 and H6. The acetylated N-terminal 
end is positioned in the interior of a triangle formed by 
residues Metl20 (H3), Pro214 (H4) and Trp255 (H5). A 
low-occupancy hydrogen bond is formed between the Ace 
CO and Tyrl21 NH. The Phe side-chain participates in a 
predominantly hydrophobic cluster, comprised of residues 
Leull7 (H3) and the non-polar moiety of Arg206 (H5) on 
one side, and Tyrl21 (H3), Alal64 (H4), Leul67 (H4), 
Leu209 (H5), Trp213 (H5), on the opposite side. 

The Orn2 mainchain carbonyl group forms a hydrogen 
bond with the NHl/2 atoms of Arg206. This hydrogen 
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Figure 3 The structure of binding mode 22. Panel A: The C5aR backbone is sliown as a thin tan tube. The first 7 residues are omitted as in 
the Nil<iforovich model [25]; the Ca atom of the 8* residue is depicted as a vdW sphere. Protein segments in contact with PMX53 are indicated 
in the following colors: 113-121 (H3) in green, 164-172 (H4) in orange, 188-199 (EC2) in black, 202-214 (H5) in yellow, 255-265 (H6) in purple, 
278-282 (EC3) in blue, and 285-286 (H7) in pink. Panels B and C present a close-up view of the binding mode and interactions. Panel C is rotated by 
180° around a vertical axis, with respect to panel B. Domains H6 and EC3 are omitted from panel B, and domains H3 and H4 are omitted from panel 
C, for clarity. H4 and H5 are additionally denoted in panels B and C to assist the reader In panels B and C, PMX53 atoms are shown in thick licorice; 
selected C5aR mainchain and side-chain heavy atoms interacting with PMX53 are shown in thin licorice. All atoms are colored by atom type. 
Hydrogen-bond interactions are shown as black dashed lines. 



bond is fully retained in the simulations of mode 22 
(100% occupancy), and is the strongest polar interaction 
between the N-terminal moiety of PMX53 (residues 1-3) 
and C5aR. An analogous interaction has been proposed to 
form in the C5aR:C5a complex, between the Arg206 side- 
chain and the C-terminal carboxylate of C5a [58]; in the 
Nikiforovich model of the C5aR:C5a complex [25] the 
Arg206 side-chain interacts with the C-terminal carboxyl- 
ate of C5a. On the other hand, the mutation Arg206Ala 
does not abrogate antagonism in linear and cyclic hexa- 
peptides [5]. It is plausible that the binding mode of these 
antagonists is somewhat different in the absence of the 
Arg206 side-chain. 

Residue Pro3 is buried between helices H5 and H6, 
interacting on the H5 side with Ala203, the non-polar 
moiety of Arg206 and Leu207, and on the H6 side with 
residues Tyr258, Met265, and to a lesser extent Gly262. 

The dCha4 carbonyl forms a medium occupancy (36%) 
hydrogen bond with the Arg206 side-chain (H5). The 
dCha4 side-chain ring participates in a cluster of "T- 
shape" interactions with residues Trp5 (PMX53), Phel72 
(H4) and His 194 (EC2). It forms additional nonpolar in- 
teractions with the non-polar part of Glul99 (EC2), and 
residues Val202, Ala203 (H5) and Met 265 (H6). 

The Trp5 ring is located at the membrane interface, 
near EC2 and helices H3 and H6. Earlier work suggested 



that this residue binds in the vicinity of Ilell6 [5]. The 
substitution Ilell6Ala correlated with the emergence of 
weak agonist activity in a PMX53 derivative with ben- 
zothiazolalanine at position 5 [5], and converted a linear 
peptide mimic firom antagonist to agonist [58]. Based on 
these results and the observation that PMX53 remains a full 
antagonist of the mutant C5aR Ilell6Ala, Higginbottom 
et al. argued that the activation of C5aR might require 
interactions between the side-chain at position 5 and res- 
idues located more deeply in the binding pocket that 
Ilell6 [5]. In their modeling of the C5aR:PMX53 com- 
plex, the proximity between Trp5 and Ilell6 was intro- 
duced a priori as a geometric restraint. The Trp position 
in binding mode 22 is consistent with these observa- 
tions and results without any restraints; Trp5 is pre- 
dicted to bind in a position less buried than Ilell6 (in 
agreement with the PMX53 antagonist activity), and in its 
vicinity (the average Ilell6 CD - Trp5 O distance is 5.4 ± 
0.6 A in the MD simulation of complex 22). Its side-chain 
forms a hydrogen bond with the Hisl94 side-chain 
and numerous nonpolar contacts with residues Cysl88, 
Glyl89, Vall90, Aspl91, Glul99. On the other side, it in- 
teracts with the side-chain of Met265. Furthermore, it 
participates in a cluster of "T-shape" interactions involv- 
ing residues Phel72 (H4), His 194 (EC2) and dCha4. Trp5 
Ca and Cp atoms interact with the non-polar moiety of 
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Figure 4 MM-GBSA interaction free-energies (l<cal/mol) for all interacting protein-ligand residue pairs in binding mode 22. Panels A 

and B correspond, respectively, to energies of non-polar and polar interactions, as defined in Eq. (3). 



ProllS and to a lesser extent with the side-chain of resi- 
due lie 11 6. 

The Arg6 mainchain CO group forms a high-occupancy 
(98%) hydrogen bond with the Tyr258 (H6) side-chain; 
the side-chain forms two persistent salt bridges with the 
Asp282 (EC3) and Aspl91 (EC2) carboxylates, a frequent 
(56% occupancy) hydrogen bond with Leu 278 O (ECS) 
and a less frequent (18% occupancy) hydrogen bond with 
Met265 S (H6). These interactions constitute the strongest 
polar interaction free energy components of mode 22 
(Figures 4 and 5). The Arg6-Asp282 interaction is in 
agreement with earlier studies, which have shown that 
Asp282 interacts with the Arg74 residue of C5a, and the 
C-terminal Arg of peptide analogs [5,59,60]. Furthermore, 
Arg6 makes numerous non-polar contacts with the non- 
polar moieties of Tyr258, Met265, Leu278, Leu282 of 
ECS, as well as nearby residues Cys287 and Val286 of H7. 
The mainchain Ca atom is placed in proximity to side- 
chain groups of ProllS and Ilell6 of HS. 



Comparison of binding modes 6, 3, 30 and 22 

The next two best binding modes were 6 and S, with cor- 
responding binding free energies of -166 kcal/mol (~9- 
10 kcal/mol larger relative to complex 22; see Additional 
file 1: Figure SI). 

Figure 6A displays a superposition of representative 
conformations from the MD runs of modes 6 and 22. 
The ligand mainchains are approximately at the same 
distance from the membrane/water interface. They are 
somewhat displaced relative to each other in a direction 
parallel to the interface (xy plane), and are rotated by ca. 
60° around the perpendicular axis Oz, and to a smaller ex- 
tent around axes vertical to planes xz and yz. The RMSD 
between the ligand conformations with the lowest binding 
free energies for modes 6 and 22 is ~ 4.5 A (ligand atoms 
N, Ca, and C). 

In the simulations of binding mode 22, PMX5S forms 
two stable intermolecular salt bridges (Arg6-Asp282 and 
Arg6-Aspl91). In the simulations of mode 6, Arg6 makes 
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Figure 5 Probability (%) maps of C5aR:PMX53 intermolecular contacts and hydrogen bonds, computed from the IVIO simulation 
trajectory of binding mode 22. Panels A, B, C and D, show, respectively, protein side-chain - ligand side-chain, protein nnainchain - ligand 
mainchain, protein side-chain - ligand mainchain, and protein mainchain - ligand side-chain contacts. A contact was considered present in a 
trajectory snapshot, if the distance betweeen the geometric center of the corresponding moieties was smaller than 6.5 A. Panel E displays a 
probability map (%) of C5aR:PMX53 hydrogen bonds. A hydrogen bond among two heavy atoms was considered present if the donor-acceptor 
distance was smaller than 3.5 A and the D-H - A angle was larger than 90°. 



Table 1 Residue pairs forming strong Intermolecular 
interactions in binding mode 22 

PMX53 C5aR^ 

Ac-Phel Tyrl21, IVlet120, Arg206, Leu117, Leu209, Trp255, Gly210, 

Pro214, Ala164,Trp213, Leul67 

Om2 Arg206, Leun7, Tyr258, Ile116, Met120, Pro113 

Pro3 Arg206, Tyr258, Ala203, Leu207, Met265 

dCha4 Arg206, Phe172, Glu199, Met265, Ala203, Val202, His194 

Trp5 Vall90, Hisl94, Glyl89, CyslBB, Asp191, Glu199, Pro113, 

IVlet265, Phel72 

Arg6 Asp282, Aspl91, Met265, Leu278, Leu281, Tyr258, llel 16 

^Intermolecular energies were computed with Eq. (3) and were averaged over 
the simulation trajectory; only pairs with total energies of at least 1.5 kcal/mol 
are tabulated. In each line, the C5aR residues are listed in descending order of 
total interaction energy with the corresponding PMX53 residue. Italics denote 
C5aR residues engaged in strong polar interactions. 



the first bridge in segment 0-9.5 ns, and the latter in the 
second half of the simulation (not shown). Numerous new 
interactions are observed: The Trp5 side-chain forms inter- 
actions with the side-chains of residues Cysl09, Leull2, 
Proll3. PMX53 segment 1-3 interacts with helix H6 in 
mode 6, and helices H3 and H5 in mode 22. 

Figure 6B displays a superposition of representative 
PMX53 conformations from binding modes 3 and 22; 
the RMSD between the conformations with the lowest bind- 
ing free energies is ~ 4.4 A (ligand atoms N, Ca, and C). 
The peptide binds at a similar distance from the mem- 
brane interface as in modes 6 and 22, and at a similar 
orientation (with its N-terminal end pointing toward the 
membrane center). Compared to mode 22, it is rotated 
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Figure 6 Comparison of PiVIX53 structures from the MD simulations of modes 6 (gray; panel A), 3 (gray; panel B), and 30 (gray; panel 
C) superposed against the lowest binding free energy conformation of complex 22 (red; all panels). The conformations shown have the 
lowest binding free energy among all snapshots of the corresponding trajectory (in the MM-PBSA approximation). The structures are aligned 
using the mainchain transmembrane atoms of C5aR. All ligand hydrogen atoms are omitted for clarity. The C5aR segments interacting with the 
ligand are shown in different colors, in tube representation; the coloring scheme is the same as in Figure 3. 



by -90° around the z-axis. As a result, the PMX53 side- 
chains point toward opposite directions and the moiety 
1-5 interacts mainly with helices H5 and H6. Arg6 forms 
polar and non-polar interactions with domains E23, H4 
and EC2, and a stable salt bridge with Glul99. Interest- 
ingly, Higginbottom et al. have argued that Glul99 might 
form an interaction with the terminal Arg side-chain of 
C5 and peptide mimics, and introduced a Glu 199 - Arg6 
salt bridge as a restraint in a C5aR - PMX53 docking 
study described in Ref [5] (Figure five(C) of that work) 
[5]. In our mode 3 this interaction is predicted without 
any guiding restraint. We note that the extracellular loop 
EC2 of that Higginbottom model had a more open con- 
formation and the Cysl09 - Cysl88 disulfide bridge was 
absent; however, the disulfide bridge is a covalent figure 
and should still be present in the open-loop conformation. 

Figure 6C displays a superposition of representative 
PMX53 conformations from binding modes 30 and 22. 
Unlike modes 22, 3 and 6, here the N-terminal end of 
the ligand is oriented toward the membrane interface 
and the N-terminal Ac group is positioned near the EC2 
loop. In this orientation, the ligand N-terminal Ac group 
is more amenable to substitutions, in agreement with 
the experimental observation that the replacement of 
the Ac group by a range of groups with variable size, 
hydrophobicity and hydrophilicity does not affect the af- 
finity or antagonist activity of PMX53 against C5aR [61]. 
Despite the change in orientation, the ligand makes the 
same two key intermolecular electrostatic interactions 
observed in mode 22 (see Additional file 1: Figure S2): 
The Arg6 side-chain forms persistent salt-bridges with 
Asp282 and Aspl91 throughout the simulation, as was 
observed in mode 22. The Orn2 CO group makes a stable 
hydrogen bond with the Arg206 side-chain in the segment 
9.5 ns - 20 ns. The binding free energy of mode 30 
is -155 kcal/mol in the first 9.5 ns and -165 kcal/mol in 
the last 7.5 ns of the simulation. Interestingly, in mode 30 



the Trp5 side-chain is positioned between the side-chains 
of Ilell6 and Val286; this pair has been suggested to form 
an activation switch for C5aR [58]. The average distances 
Trp5 CH2 - Ilell6 CG2 and Trp5 CH2 - Val286 CG are, 
respectively, 4.6 ± 1.2 A, and 5.0 ± 0.6 A. 

Comparison of binding mode 22 and ttie C5a:C5aR complex 
(Nildforovicti model) 

A detailed comparison of the promising binding modes for 
the C5aR:PMX53 complex, predicted in this work, and the 
Nikiforovich model for the C5aR:C5a complex [25] can 
help interpret the antagonist activity of PMX53. Figure 7 
displays the binding conformations of C5a and PMX53 li- 
gands, obtained by alignment of the intra-membrane C5aR 
mainchain heavy atoms in the C5aR:C5a complex (Nikiforo- 
vich model) and C5aR:PMX53 complex (mode 22) with the 
program SuperPose [62]. Additional file 1: Figure S3A-D 
display maps of side-chain and mainchain contacts in 
the C5a:C5aR Nikiforovich model; [25] intermolecular 
hydrogen bonds in the C5a:C5aR complex are presented 
in Additional file 1: Figure S3E. 

According to the superposition presented in Figure 7, 
PMX53 blocks the entry of the C5a C-terminal into the 
transmembrane domain of C5aR, by occupying the pos- 
ition in which the C-terminal domain of C5a binds and 
promotes signaling. PMX53 does not directly act as a 
mimetic of the C5a C-terminal, as PMX53 is cyclic and 
the C-terminal domain of C5a is linear. Upon superpos- 
ition of the PMX53 and C5a in complex with C5aR, the 
N-terminal end of PMX53 coincides with the C-terminal 
end of C5a, while the 71-73 C-terminal domain of C5a 
is placed between the PMX53 3-5 backbone moiety on 
one site and the Orn2 side-chain moiety on the other 
site. Also, the backbone of C5a residues Asp369 and 
Met370 structurally coincide, upon superposition, with 
the side-chain of PMX53 residue Trp5. The C5aR resi- 
due Aspl91 makes a salt bridge with PMX53 residue 
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Figure 7 The C5aR:PMX53 complex (conformation of lowest 
binding free energy in the simulation of binding mode 22) is 
superposed against the Nikiforovich model [25] for the C5aR: 
C5a complex. The C5aR mainchain (mode 22) is indicated by a thin 
tan tube, and the Ca atom of the first N-terminal end residue of the 
simulation system (8* residue of C5aR) by a vdW sphere. Interacting 
protein domains are colored in different colors; the coloring scheme 
is as in Figure 3. C5a is displayed as a red tube and the Co atom of the 
1^' residue is indicated by a vdW sphere. PMX53 is shown in licorice, 
colored by atom type. Hydrogen atoms are omitted for clarity. 



Arg6 in mode 22, and with C5a residue Lys68 in the 
C5aR:C5a complex. Furthermore, the first four PMX53 res- 
idues interact strongly with H5 residues Ala203, Arg206, 
Leu207, Leu209; these interactions could compete with 
interactions within the C5a 69-74 residue moiety, and 
particularly with the salt bridge between Arg206 and the 
C-terminal carboxyl group of Arg74. Last but not least, 
H6 residues Gly262 and Met265 interact with C5a seg- 
ment 67-70 and Arg74 in the Nikiforovich model; these 
interactions may also be blocked by interactions of C5aR 
with Orn2, Pro3, dCha4 and Trp5 of PMX53. 

Conclusion 

In the present work we construct structural models for 
the complex of membrane-embedded C5aR and its an- 
tagonist peptide PMX53 via a computational framework 
that combines conformational sampling for both recep- 
tor and ligand, docking and filtering of conformations by 
structural and energetic criteria. A large body of experi- 
mental results suggests that the C5aR binding site of the 
antagonist cyclic hexapeptide PMX53 is at or near the 
TM binding site of the C5aR agonist peptide C5a [5]. 
The key PMX53 residue Trp5, an important determinant 
of antagonism, is likely to be positioned near Ilell6, a 
residue implicated in interactions with the PMX family 
of peptides and possibly the activation of C5aR [58]. Fur- 
thermore, Arg6 has been shown to interact with Asp282 



[5]. The most promising binding mode (22) and several 
other modes of low binding free energy reproduce both 
interactions. An additional interaction is observed in 
several high-affinity modes, between the mainchain CO 
group of Orn2 and the Arg206 side-chain. An analogous 
interaction has been proposed to form in the C5aR:C5a 
complex, between the Arg206 side-chain and the C- 
terminal carboxylate of C5a [58] and is present in the 
Nikiforovich model of the C5aR:C5a complex [25]. 

Notably, the aforementioned interactions are predicted 
without imposing any a priori geometric constraints in 
the initial conformations of the complexes. This success- 
ful outcome lies in the use of a multi-step computational 
framework, which included the generation of representa- 
tive receptor and ligand conformations by long MD sim- 
ulations with high quality implicit-membrane models, 
the docking of a large combination of receptor and ligand 
conformations, the comprehensive filtering of several hun- 
dred thousand conformations for the complex by struc- 
tural and free energy criteria, and the re-evaluation of the 
most promising binding modes by additional MD simula- 
tions and binding free energy calculations. 2. Membrane 
effects were introduced by high-quality implicit-membrane 
models [29-31], which enabled both the rapid generation 
of representative structures and the estimation of their 
binding affinities. Therefore, the computational framework 
used here is capable of addressing the challenges in gener- 
ating reliable structural models for GPCR protein-ligand 
complexes of unknown structure. While, the MM-PBSA 
and MM-GBSA methods yield large binding free energy 
values [53,54,63], their use in the specific computational 
protocol of this study proved useful, with regard to identi- 
fying the lowest binding free energy, and thus, the most 
promising PMX53:C5aR binding mode. The large affinity 
free energies are partly attributed to the omission of the 
protein, ligand and complex configuration entropy contri- 
butions to binding; due to energy-entropy compensation, 
when these terms are included in the calculation, they are 
expected to yield significantly smaller total free energies 
[64]. These entropic terms are associated with large uncer- 
tainties and are expected to cancel out to a large extent in 
the relative affinities of different binding modes. Therefore, 
their omission is not expected to be important at the level 
of accuracy of the present calculation, which identifies 
most promising binding modes. In a similar fashion, the 
membrane MM-PBSA and/or MM-GBSA approximations 
were used to identify the most promising binding modes 
of a dual tropic HIV-1 gpl20 V3 loop in complex with 
CXCR4 and CCR5, CXCL12 (SDF-la) in complex with 
CXCR4 and CCL5 (RANTES) in complex with CCR5; the 
computationally derived structures were in exceptional 
agreement with experimental findings [65-68]. 

The structure of membrane-embedded C5aR:PMX53 
complex and its dynamic features, presented here, will 
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serve as a template for biopharmaceutical discovery of 
peptide, peptidomimetic, and organic compound antago- 
nists of C5aR for targeting of complement-mediated auto- 
immune and inflammatory diseases. Such antagonists may 
have superior ADMET (absorption, distribution, metab- 
olism, excretion, and toxicity) properties compared to 
PMX53 and therefore be more amenable to clinical devel- 
opment. Also, the methodology developed for the gener- 
ation of the structure of membrane-embedded C5aR: 
PMX53 prepares the grounds for the generation of similar 
models for the homologous receptors C3aR and C5L2 
with their respective endogenous ligands, as well as ago- 
nists and antagonists. C5aR, C3aR, and C5L2 have com- 
mon as well as distinct structural features and biological 
functions. Comparative studies at atomic resolution will 
shed light on the underlying structural, physicochemical, 
and dynamic properties that mediate the similarity and 
variability in their structures, dynamics, binding proper- 
ties, and biological functions. 

Although the database of solved GPCR structures is 
slowly but steadily increasing, and there are more struc- 
tural templates available today for computational modeling 
than previously, still there is no structural information for 
the vast majority of GPCRs; and structural information on 
GPCR-ligand binding at atomic resolution is rare. The 
computational framework presented here can be of wide 
use for the development of GPCR-ligand structural models 
in membrane environments. Such models will be useful in 
providing the structural basis for mechanistic studies of 
the interactions between GPCRs and ligands, and their ef- 
fects in intra-membrane and intra-cellular dynamics that 
drive the selection of intra-cellular activation pathways. In 
addition, such models will be useful in knowledge-based 
biopharmaceutical discovery. 
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